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^ Abstract We present a discretization of the linear advection of differential forms on 

•/^ bounded domains. The framework established in [4] is extended to incorporate the 

Lie derivative, X, by means of Cartan's homotopy formula. The method is based 

on a physics-compatible discretization with spectral accuracy . It will be shown that 

the derived scheme has spectral convergence with local mass conservation. Artificial 

dispersion depends on the order of time integration. 



1 INTRODUCTION 



Consider the classical advection problem for a scalar function in conservation form, 



^ do 

\0 Jl+V-(vp) = 0, (1) 



where v is a prescribed uniformly Lipschitz continuous vector field and p the ad- 
vected scalar function. The method presented in this work is based on the approxi- 
mation of the differential operators with the focus on the spatial discretization and 



-y-s a time-stepping scheme that distinguishes between quantities evaluated at time in- 

^— I stants and quantities evaluated over time intervals. 

L| The mimetic framework, presented in [4], showed that using the differential ge- 

ometric approach for the representation of physical laws clarifies the underlying 
structures. One clearly identifies to what kind of geometrical object a certain phys- 
ical quantity is associated and this determines how its discretization must be done 
(e.g.: evaluation at points, integration over lines, surfaces or volumes). Additionally, 
a well defined, metric free, representation of differential operators is obtained, to- 
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gether with their metric dependent Hilbert adjoints. For these reasons, the authors 
followed this approach for the advection equation. It is known, see [1, pp. 317], that 
(1) is a particular case of the generalized advection equation which can be written 
in terms of differential geometry as, 

^.Xva^^^=0. (2) 

at 

The advection operator, Xv, is the Lie derivative for the prescribed velocity field 
V and the advected quantity is given by the ^-differential form a^'^\ Depending on 
the index k the quantity a*^*' can represent scalar, vector and higher dimensional 
quantities. 



2 DIFFERENTIAL GEOMETRY 

In this section a brief introduction to differential geometry is given. For a more 
detailed introduction the reader is directed to [1]. Given an n-dimensional smooth 
orientable manifold Q it is possible to define in each point a tangent vector space 
E of dimension n. The space of smooth vector fields on a manifold is the space, F, 
of smooth assignments of elements of E to each point of the manifold. We denote 
by A*^, k an integer < A: < n, the space of differentiable A:-forms, i.e. all smooth 
^-linear, antisymmetric maps cJ-'^'' : £ x ■ ■ ■ x £ — > R, at every point of the manifold. 
We recall the wedge product A : A* x A' — > A*^"""' for ^ + / < n with the property that 
Q,W /\^(0 - (^-i^kipii) ^ ^(i) -pjjg jj^j^gj- product (•,•) on E induces at each point of 
the manifold an inner product (■,) on A'. In turn, this can be extended to a local 
inner product on A*^, [8, pp. 149]. The local inner product gives rise to a unique 
metric operator, Hodge-*, • : A*" -> A""*^, defined by a**' A •yS^'"' = (a**^', /?**') a)'-"\ 
where cJ-"^ = •I is the standard volume form. By integration, one can define an inner 
product on Q as (■,-)l2 '■- Jni''')'^'^"^- The exterior derivative d : A*^ — > A*^^^ satisfies 
the following rule, d(a'(*' A/?®) = da* A^''' + (-l)*^^^ a d/3<" and by definition 
da^"^ - 0. The flat operator, b, is a mapping b : T i-> A'. 

The Lie derivative along a tangent vector field, v, is denoted by Xv and represents 
the advection operator in differential geometry. It is a mapping Xv : A*^ i-> A^. From 
Cartan's homotopy formula the Lie derivative can be written as 

XvO'^'''():=diva'('^' + tvda*, 

where the interior product of a tangent vector field, v, with a A:-form, a^''\ is a 
mapping lyO'''*^ : A*^ — > A*^"' given by: 

tya^''\X2,---,Xi):=cK<*'(v,X2,---,Xi), VX; e T and i^a^'^^^O, Vv e T . 
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The interior product is the adjoint of the wedge product, made explicit by: 

where v*" = v*'' e A^ and a^*^' e A^. 

h 



The relevance of this adjoint relation between the interior product and the wedge 
product lies in the fact that it shows how a physical quantity represented by an 
interior product with a vector field can be represented by its dual differential 1-form. 

For a volume form p^"' the Lie derivative is simply XvP*^"' = divP*"^ and for a 
0-form, Xva'^'^ivda*"^. 



3 MIMETIC DISCRETIZATION 

In this section a brief introduction to the discretization of physical quantities and to 
the discretization of the exterior derivative is presented. For a more detailed presen- 
tation the reader is directed to [2, 4, 6] 

Consider a three dimensional domain Q and an associated grid consisting of 
a collection of points, T(0),,, line segments connecting the points, T(i),,, surfaces 
bounded by these line segments, T(2),i, and volumes bounded by these surfaces, T(3) ,-. 

Let A^ be the space of smooth diflFerentiable A:-forms. Additionally, let 
the finite dimensional space of differentiable forms be defined as A^ = 
span({e. )), i- !,■■■ ,dim(A^), where e. € A*^ are basis A:-forms. Under these 
conditions it is possible, see [4, 6], to define a projection operator tt/, which 
projects elements of A* onto elements of A^ which satisfies: 



It is possible to write: 



TThd = dnh ■ (4) 



(k) (k) V (*:) 



where 



ik) V ( 

r a^*^) and f ej 



afi= I a-"' and | ef^Sij, fe = 0,l,---,n 



A set of basis functions yielding a projection operator iih that satisfies (4) can be 
constructed using piecewise polynomial expansions on the quadrilateral elements 
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using tensor products. Thus, it suffices to derive the basis forms in one dimension 
on a reference interval and generahze them in n dimensions. 

In one dimension take a 0-form, a^''' e A^(Qref\ where Q^gf :- [-1, 1]. Define 
on Qref a cell complex D of order p consisting of (/? + 1) nodes T(0)_, = ^, with 
i -Q,--- ,p, where - 1 <^i) < ■■■ <^i < ■■■^p < 1 are the Gauss-Lobatto quadrature 
nodes, and p edges, T(i),- - \^i-\,^i\ with / = I,--- ,p. The projection operator jih 
reads: 

p 
7rAa(">(^) = ^«,ef©, (5) 

i=0 

where e (^) = /,- (f ) are the p*^ order Lagrange polynomials and a,- = cfii^d- Simi- 
larly in one dimension for the projection of 1-forms Gerritsma [3] and Robidoux [7] 
derived 1 -form polynomials czWeA edge polynomials, e. e AHQref), 

i-i ,, 
ef'© = e,©d^, with eiiO^-Y-^. (6) 



Note that in this way we have: 

Moreover, the exterior derivative of the basis 0-forms is given by: 






(7) 






K k=0 ^ 



= 6P-e||>, i=l,---,p-l. (8) 



In this way, the exterior derivative of a discrete 0-form can be written as: 

d«f = dX«,^"'= t ^^a,^. (9) 

,=0 i=0,./=l 

where, E..' is the incidence matrix containing only the values, 0, 1 and -1, see 
[2, 4, 6] for more details. This idea can be extended to higher dimensions, giving 
rise to ^-incidence matrices, E . ^ ' , which represent the discrete exterior derivative 
on discrete A:-forms, see [2, 6]. 



4 Mimetic spectral advection: an application to ID advection 

In this section we want to illustrate how to discretize the advection equation. Take 
the Lie advection of a 1 -form. 
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at 



at ~ "^ 

(10) 



Here p^^^ is the advected quantity, say mass density, and ^^^^ represents the instanta- 
neous fluxes of the advected quantity under the vector field v, which are discretized 
in space as, 

^-t^^i^^ and p('> = |]p,(0.r' . (11) 

1=0 1=1 

For the sake of clarity in the method presentation we first introduce the time treat- 
ment, then the interior product discretization and finally their combination for a 
numerical solution of the advection problem. 



4.1 Time integration 

The time integrator used for solving the time evolution part of the advection equation 
is the canonical mimetic one, an arbitrary order symplectic operator derived in [5], 
which is connected to canonical Gauss collocation integrators. Take an ordinary 
differential equation of the unknown function y(f): 

— =/i(y,0, fe/cR. (12) 

df 

Discretizing y{t) as yh = Tjl^Q^hit), one gets: 

*:=0 k=l 

where the superscript k denotes the time level. 

The approximated solution, yh{t), is a polynomial of order p determined by 
means of {p + 1) degrees of freedom such as its values at the Gauss-Lobatto nodes, 
red dots in Figure 1. On the other hand, -g^ is a polynomial of order {p - 1) defined 
by only p degrees of freedom. One can set these degrees of freedom to be the values 
of the derivative in one point inside each of the p intervals [f*^,?*^"*"']. A choice that 
results in a symplectic integrator of order 2p is to select these points as the Gauss 
nodes of order (/? - 1), the blue nodes of Figure 1. Notice that along the trajectory 
these nodes will not show the usual Gauss-Lobatto and Gauss distribution patterns, 
since in general the velocity field is not constant. In this way the discrete integrator 
becomes: 
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with V the p nodes of a Gauss quadrature formula. The fact that the instants in 
time, t'^, where the y^ are defined alternate with the instants in time, P, where the 
hi are evaluated (see Figure 1), corresponds to a staggering in time. This staggering 
also appears in leap-frog methods and in the implicit midpoint rule, for instance. 





-- ^ ^ L 


f f / / . v^ 

Uj'- 


-- ^ ^ L 


f f / / . 


-^ ^ ^ L 

^ ^ / f 
. . / f f 



rig. 1 Geometric interpretation of the solution of (12) as given by (13): (^^'"'(f))- In red the Gauss- 
Lobatto nodes where the trajectory is discretized. In blue, the Gauss nodes where its derivative is 
discretized. The flow field, represented by arrows, is tangent to the curve at the Gauss nodes. That 
is, the derivative of the approximate trajectory is exactly equal to the flow field at the Gauss nodes. 



The first equation in (10) using the discretization (11) and (9) can be written as: 



df 



(1) 



^-Y^m^ 



dpi(t) 
dt 



^-Y,^\,(t). (14) 



This equation has a similar form as (12), but now as a system of equations, therefore 
one can apply the mimetic integrator, yielding: 



X(pf--p?)..(?^)=-2Ei;'V 



(15) 



Recall that p^ is the discrete degree of freedom of the advected quantity at the t'^ 
instants of time associated to Gauss-Lobatto nodes and ^ J' is the discrete degree of 
freedom of the fluxes of the advected quantity at the ?' instants of time associated to 
the Gauss nodes, just as stated for the systems of ordinary differential equations. 
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4.2 Interior product 

The discretization of the interior product is done using (3), in the following way: 



Definition 1 (Discrete interior product). 

terior product ty,/, : A^ — > AJJ is such that: 


In one 


dimension, 


the discrete in- 


(...»<".<;»'),. 


= K^V^A 


(0)\ 


Sef^A 


i (16) 


where v'' = v<''eA^ and a^^' € 


^\- 









In this way one satisfies the duality pairing between the interior product and the 
wedge product in the discrete setting. 

Partitioning the domain £? in a spectral element cell complex one can apply the 
discretization of the interior product in each spectral element, obtaining: 

Yj>^{^^^''^^%-Y^^^^^^ ' ^^'-^- (17) 



4.3 Putting things together: advection 

The complete discrete systems becomes: 

(18) 



5 Numerical results 

This approach was applied to the two dimensional solution of an advected sine wave 
and a sine bell in a constant velocity field v = e^: p^^\x,y) = sin(;7rx) sin(7r3')dxdy 
(sine wave) and p^^\x,y) - sin(27rx) sin(27ry)dxd}' if (x,y) e [0,0.5] x [0,0.5] and 
p^^\x,y) = in (x,y) e R^\[0,0.5] x [0,0.5] (sine bell), on a domain with periodic 
boundary conditions. 

In Figure 2 the error in time of the numerical solution of ( 1 0) for a mesh of 4 x 4 
elements with a zff = 0.1.? and various polynomial orders in space, p, and time, pt, 
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is presented. The initial error, due to the discretization, is conserved, as long as the 
time integration is sufficiently accurate. 
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Fig. 2 Error in time of the numerical solution of (10) with v = e.j and 4x4 elements of order p = 3, 
p = 6, p = 10 and p = 12 (from left to right and top to bottom) and At = O.ls, for the sine wave 
p^^\x,y) = sin(;rx) sin(;ry)djcd>'. As shown, the error in the solution increases with time due to the 
inaccuracy of time integration. When time integration is accurate enough the error in the initial 
state is preserved. Here p, denotes the polynomial degree in time. 



In Figure 3, the h- and p-convergence plots are shown for different values of the 
order of the time integration scheme, pt, and At = 0.1.?. It is possible to see that the 
method presents algebraic /j-convergence rates of order (p+l) as long as the time 
integration error does not dominate the spatial one. The method shows a spectral 
p-convergence as soon as the time integration is accurate enough. 

In Figure 4 the error on the velocity is presented as a function of the advected 
sine wave frequency. This figure shows that the numerical method introduces an 
artificial dispersion if the time scheme is not accurate enough. 

Another fundamental aspect is the conservation of the advected quantity. Figure 5 
shows the mass error in time, that is: Lpj - LpL • The error goes from the zero 
machine in the first 10^ time steps while thereafter it steadily increases. Notice that 
even after 2 x 10"^ time steps the error is still below 10"^^. 

In Figure 6 one can see the advection of a sine wave of frequency w = ;7r in a Rud- 
man vortex for 100 time steps after which the direction of the flow is reversed and 
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Fig. 3 Left: h convergence in space for the advection of a sine wave with At = O.ls. Right: p 
convergence in space for the advection of a sine wave, 4x4 elements and At = Q.ls. 




Fig. 4 Error in velocity as a function of the frequency of the advected sine wave: numerical dis- 
persion, p = 10, At = 0.1.V and ;i = 4x4 elements. 




Fig. 5 Sum of the local errors of the advected 2-form for a sine bell in a velocity field v = e.,, with 
50x50 elements of order p = (blue) and 4x4 elements of order p = 9 (red),zlf = 0.01 and pt = 2. 



the calculation is continued for another 100 time steps, with 4x4 curved elements of 
order p -9,At = Q.ls and time integration of order /?, = 2, on a distorted mesh. The 
mimetic advection enables one to recover the initial solution, thus demonstrating 
that the integration method is reversible. 
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Fig. 6 From left to right and from top to bottom: advection of a sine wave of frequency w = ;t on a 
Rudman vortex with 4x4 curved elements of order p = 9, At = 0.ls and time integration of order 
p, = 2, on a distorted mesh. At time f = 10, the flow field is reversed. At times t = 8.0.V and t = 12.0.V 
the mesh is visible in the solution. See http : / fwwv . youtube . com/watch?v=Qmo Jyqtk9YA for 
animation. 
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